APS/123-QED 



O 

> 
O 



43 



van't Hoff-Arrhenius Analysis of Mesoscopic and Macroscopic Dynamics of Simple 
Biochemical Systems: Stochastic vs. Nonlinear Bistabilities 
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Multistability of mesoscopic, driven biochemical reaction systems has implications to a wide range 
of cellular processes. Using several simple models, we show that one class of bistable chemical sys- 
tems has a deterministic counterpart in the nonlinear dynamics based on the Law of Mass Action, 
while another class, widely known as noise-induced stochastic bistability, does not. Observing the 
system's volume (V) playing a similar role as the inverse temperature (f3) in classical rate theory, 
an van't Hoff-Arrhenius like analysis is introduced. In one-dimensional systems, a transition rate 
between two states, represented in terms of a barrier in the landscape for the dynamics &(x, V), 
k oc exp{ — VA$* (V)}, can be understood from a decomposition A&(V) ~ Ac/>J + A<f>\/V. Non- 
linear bistability means A<^jjj > while stochastic bistability has AcAg < but A(f>\ > 0. Stochastic 
bistabilities can be viewed as remants (or "ghosts") of nonlinear bifurcations or extinction phe- 
nomenon, and A0q and AcA* as "enthalpic" and "entropic" barriers to a transition. 

PACS numbers: 87.10.-e;64.70.qd;02.50.Ey 
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Biochemical reaction dynamics in a small volume on 
the order of a cell is stochastic. For a given spatially ho- 
mogeneous biochemical kinetic mechanism, be it a gene 
regulatory network or an intracellular signaling path- 
way, the stochastic trajectories of the chemical compo- 
sitions of a mesoscopic, nonlinear reaction system can 
be computationally modelled via the Gillespie algorithm, 
and its probability distribution follows the chemical mas- 
ter equation (CME) first studied by Delbriick Q. The 
Dclbruck-Gillcspie process (DGP) is a multi-dimensional 
birth-and-death process [2[, with the system's volume, 
V, as a key parameter. In the limit of infinite large V, 
a macroscopic nonlinear dynamical system emerges §■ 
This is precisely the system of ordinary differential equa- 
tions (ODEs) following the classic Law of Mass Action 
(LMA) for chemical kinetics. 

Bistability in terms of two stable fixed points is one of 
the salient features of many nonlinear chemical reaction 
systems The stationary distribution of a DGP that 
corresponds to a macroscopic bistable nonlinear chem- 
ical reaction system is bimodal 5]. Recently, it also 
has been discovered that certain nonlinear reaction sys- 
tem with small V can exhibit bimodal stationary distri- 
bution which has no macroscopic bistable counterpart. 
This phenomenon has been called noise-induced bistabil- 
ity and stochastic bifurcation 0, 0] ■ Recent experiments 
in synthetic biological systems have partially confirmed 
the theoretical insights Q. 

One of us has proposed the notion of nonlinear bista- 
bility and stochastic bistability to distinguish the two dif- 
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ferent scenarios 0. In particular, it was suggested that 
there is a distinct difference in the volume dependence 
of the transition rates between states. While the transi- 
tion rate of a state decreases with increasing V for the 
nonlinear bistability, it actually increases for stochastic 
bistability. A van't Hoff-Arrhenius-like analysis with re- 
spect to system size V seems possible Q- 

For a large class of DGP, the stationary probability 
distribution with increasing V has the asymptotic ex- 
pression [T3] 



Pn(V) 



-V(t>a(x)—<)>\{x) 



(1) 



where x — n/V is the concentration(s) of the chemical 
species, and the <p's are independent of V. Forthermore, 
it can be shown that if <j)o(x) exists and is differentiable, 
then it is a Lyapunov function for the macroscopic ODE 
dynamics dx/dt = F(x) fill ]: 



dt 



< 0. 



(2) 
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Therefore, the stable fixed points of the ODE are located 
at the minima of </>o(x), and the leading term in the ex- 
ponent in Eq. (JXJ) indicates that the "barrier" between 
two stable fixed points increases with V. In the limit 
of V = oo, ergodicity breaks down and there will be no 
transitions between stable fixed points (attractors). 

One can making an analogue between V~ x and temper- 
ature T in the traditional rate theory Q. Both V^ 1 — >• 
and T — > imply a deterministic limit. In fact, Ar- 
rhenius law states that a rate constant k oc e _AG / fcfiT , 
where the activation free energy, according to van't Hoff 
analysis, has an enthalpic and an enthropy part AG* = 
AH* — TAS$. Therefore negative activation enthalpy 
leads to decreasing k with increasing temperature [l2j . 
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Comparing the van't Hoff-Arrhenius analysis with Eq. 
([1]), we can identify 0o with "enthalpy" and 0i with en- 
tropy. As we shall show below, stochastic bistability is 
associated with a negative A0J. 

The general theory — Let us consider the 1- 
dimensional birth-and-death process with birth rate 
u n iV) and death rate w n (V): 



— p n (t) = Pn-lU n -i -p n (w„ + U n ) + p n+1 W n+1 , (3) 

at 

(n > 0,u_i = wq = 0.) The unique stationary distribu- 
tion of Eq. ([3]) is 



Pn — Po 



n 

1=0 



u t {V) 
w e+1 (V) 



:pg t e -V»(n,V) j (4) 



where pff is a normalization factor, and <&(n, V) 

n-l 

V 



1 71 — 1 



e=o 



we +1 (V) 



u t {V) 



0o(x) + -^(z,n (5) 



here we have assumed n — xV. To obtain 0o(x), we let 
V — > oo while holding x constant. 

To derive Eq. ([5]) , we consider £ = V z and noting that 
functions u and w have asymptotic expansions according 
to the macroscopic LMA 



V- 1 u Vz {V)^ f jL (z) + V- 1 f jL 1 {z) + - 
V~W»(V) « A (z) + l^Aifc) + 

Then we have: 



0o (x) = 



(6) 



(7) 



We can also obtain a leading order approximation for 
ip(x, V) w </>i (x) + V _1 02(a) H ■ Note that 



o V z - 



xV-1 

E 



In 



J_ f x d_ 

ZV J d~z 



Xojz) 
Mo(z) 

A Q (z) 



(=0 



A (l/VQ 



« 7 < 



V 



V 



In 



and 



E ln 



Mo(z) 
u>*+i(V) 



dz + o(V- 1 ) , 



iV-l 



E ln 

1=0 



X (£/V) 



x fX 1 (z)+X' (z) 



Ao(«) 



«>(z) 



dz + o(l) . 



Therefore, 

01 (x) = 



Ai(z) _ MjXgA rfz + ln(/x (x)A (x)) 
o VAo(^) Mo(z)/ " 2 



Therefore, $(x) « 0o(x) + (f>t{x)/V. 

Stochastic bistability — Stochastic bistability means 
for a finite V, the $(x, V) has a minimum at x = x* 
and a maximum (or saddle point in multi-dimensional 
problems) at x*, with $(x*, > $(x*,V), but 0o(x :t ) < 
0o(x*). Noting the relation $(x, V) « 0o(x) + V^ _1 0i(x), 
this indicates that 0i(x : ' : ) — 0i(x*) = Acf)\ > 0. 

We illustrate the theory by an example. The Schlogl 
model of nonlinear chemical reactions with autocatalysis, 



Oil Pi 

A + 2X=±3X, X =± B, 

02 02 



(9) 



has recently found wide applications in cellular biochem- 
istry Q. With appropriate parameters it is well-known 
to exhibit bistability. We now consider this model out- 
side but near its bistable regime. The stochastic DGP 
has birth and death rates 



Un(V) 
W n (V) 



kin(n — 1) 
V 

fc_in(ri — l)(n — 2) 
V 2 



(10) 



k 2 n, 



in which k± — ai[A], fc_i = a 2 , fc 2 = /?i and fc_ 2 
& [-B]- Then, /«o(^) = fciz 2 + fc_ 2 , ^i(z) = —kiz, \o(z) 
k_ 1 z 3 + k 2 z, Xi(z) = -3fc_iz 2 . Then, 



0o (x) = x In 



wy)x „ /i/ 

— 2a— arctan 



x, 



01 (x) 

in which 





(11) 
(12) 



8 = ki/k 2 , v = k- 2 /k 2 , 7 = fcifc 2 /(fc-ifc_ 2 )- (13) 

Fig. [1] shows that with increasing V, the bistability dis- 
appears in the deterministic limit. Furthermore, it shows 
that the "activation enthalpy" A0j has a negative value, 
and the barrier at finite V is an "entropic" one. This is 
the origin of the stochastic bistability. 

The ghost of extinction — A canonical phosphory- 
lation-dephosphorylation signaling with positive feed- 
backs exhibits stochastic bistability Q: 

E + K + E* ^ 2E* + K, E*+P^E + P, (14) 

in which A" is a kinase that catalyzes the phosphorylation 
reaction E — ?> E* , and P is a phosphatase that catalyzes 
the dephosphorylation reaction. The K is only active, 
however, after binding an E* [TJ. The rates of the DGP 
are 



km(N 



V 



-+k-2(N—n), w r , 



k-\n{n — 1) 
V 



-k 2 n, 
(15) 
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FIG. 1: Stochastic bistability has a barrier height A&(V) 
which decreases with increading V, as revealed by the van't 
Hoff-Arrhenius analysis A$* = + V r_1 A^. (A-C): Blue 
solid line: <f>o(x) where x = n/V is concentration of X in 
Black dots: exact $(x, V); Red dashed line: </>o(x) + V~ (j>i{x) 
according to Eqs. (|llll2p with 6 = 0.25, v = 1.2, 7 = 15, and 
different Vs shown in the figures. (D) van't Hoff-Arrhenius 
plot, with filled squares for A&(V), shows a slope ~ A0* > 
and an interaction ~ A0q < 0. 



where fci = ai[K], fc_i ~ a 2 [-^]; ^2 = Pi[P]i k-2 — 
P2[P] and N is the total copy number of E and E* 
molecules together. We have ^0(2) = (fciz + fc_ 2 )(e t —2), 
/ii(z) = 0, Ao(z) = k-\z 2 + k2Z, and Ai(z) = — fc_iz, 
where ej = N/V . Therefore, 

9x + vy (9 \ 

4>o(x) = In — x + l\+x\n.x (16) 

9 \isj J 



9x + v 



\n(9x + v) + (e t — x) ln(et — x), 



h{x) = -In 



x(9x + v)(e t — x) 



9x + z/7 



(17) 



in which 9, v, 7 are again given in Eq. (I13[) . 

The stochastic bistability in © and Fig. Q] occurs 
when the system is near its deterministic bistable regime. 
Therefore, one can consider the stochastic bistability as 
a "ghost" of the saddle- node bifurcation [3]. In the 
present case, 4>o(x) in Eq. (p~6|) is a convex function with 
a single minimum on [0, et] for all parameters. Therefore, 
the deterministic dynamics of (1141) can not have bistabil- 
ity or bifurcation. It however, can have a unstable fixed 
point at x — when /?2 = 0. This is the case of "extinc- 
tion". Even though the x = is unstable to the ODE, 
the stationary distribution for the CME is pff — S n o: The 
stochastic dynamics goes extinct with probability 1. As 
a "ghost" of the extinction, this system can also exhibit 
stochastic bistability if /3 2 is nonzero but sufficiently small 
Fig. [5] shows that while <f>o(x) is convex. However, 
for finite V there is a stochastic stable state at x = 0. 



FIG. 2: Solid blue lines show 4>a(x) according to Eq. (|16|) . 
Black dots are &(x, V) following Eq. ((5|. Parameters used 
are 6 = 0.7, v = 0.001, 7 = 1000, e t = 10, with (A) V = 2 and 
(B) V = 10. x = is a stochatic stable state which disappears 
in the deterministic limit. (C) van't Hoff-Arrhenius plot with 
filled squares for A$* showing a negative A0q, as expected 
for stochastic bistability. 



System with single molecules and stochastic 
bistability — We have so far assumed that each term in 
the rates u n (V) and w n (V) corresponds to a term in 
uo(x) and \q(x). However, if a chemical reaction sys- 
tem at finite volume contains a term O(l), then in the 
limit of V~ 1 uy x (V) ~ ^0(2) + V^ 1 ^i(x) + • • • , the term 
only contributes to fj.i(x). This is in fact the so-called 
"single-molecule effect" . jl5j 

DGP with Lewis' chemical detailed balance — 
For chemical system with a single dynamical species, 
the DGP predicts that chemical equilibrium has either 
a binomial distribution (canonical ensemble) or Poisson 
distribution (grand canonical ensemble) following G.N. 
Lewis' principle of chemical detailed balance [l6j . 



The canonical ensemble has u r . 



k+(N 



n),w n = 

k-n (n < N), (f>o(x) = x\n(k_x/k + ) + (et — x) ln(et — 
x), and (pi (x) — jh\(x(e t — x)), with x < et — N/V. 
<f>oi x ) is convex with a minimum at x = k + e t /(k + + 
£;_). More importantly, V4>q(x) + <f>\(x) is convex for 
x E (^ 7 e t — ^7). Note that the asymptotic expansion 
ip(x, V) = 4>i(x) + O (V -1 ) is not uniformly valid at x — 
0, et- The existence of the ghost of extinction can not be 
determined by <$>i(x). 

The grand canonical ensemble has u n — jV, w n = kn, 
(po(x) — a;ln(a;/x*) — x, and 4>i(x) = ^lnx, where x* = 
j/k. Again, convex 4>q(x) has a minimum at x* , and 
<f>i(x) is concave. Still V4>q(x) + (/)i(x) is convex when 
(2x7,00). 

Therefore, with detailed balance in a chemical reaction 
system, the equilibrium distribution is always unimodal. 
The 7 parameter in the previous sections represents the 
energy dissipation of open chemical systems. One can 
verify that when 7 = 1, the results in the previous two 
sections are reduced to what we have here. 

Exit rate of a stochastic stable state — So far we 
have exclusively discussed the stationary distribution in 
the form of e - ^*^'^ and its relation to the stochastic 
bistability. We now establish the relation between the 
transition rate from one stable state to another and the 
function &(x, V). In a Id birth-and-death process, the 
mean time of the first arrival at n\ starting at n\, with 
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reflection at (0 < n\ < nJj), has been widely used in 
the various lattice hopping models [l7j : 



E E 



Pn(V) 



(18) 



in which p s * is the stationary distribution to Eq. ([3]). 
Let n\ = Vx\ , = Vxt, , then we found the asymptotic 
expansion formula T x j^ x * = 



r 



A(») Jo 



1 + 



V 



in which, the modified V") 
and A(x) = 



1 J , / v , , Mo(g)/Ao(aj - 1 
n Wl)+ln in(,oWAoW) 



1 



A (x) - hq{x) 
u (x) \ln X (x) - \n fi (x) 



(19) 
(20) 

(21) 



Noting A(x) playing the role of diffusion coefficient. Eq. 
(fl~9| is the same as Kramers' theory for the rate of cross- 



ing a continous energy barrier located at x* G (x^x?";) 

(22) 



me a 

17]. Applying Laplace's method leads 

2lle VmxKv)~3>{xl,V)] 

Mxi)^{xl)\4>l(xi)\ 



T 



This simple relation between transition rate and $(x, V) 
are unique for 1-d system. While the concept we de- 
veloped here will be valid in higher dimensional systems 
with saddle point, computation will be demanding [9]. 



Summary — Nonlinear biochemical reaction systems 
can have bi- or multi-stable behavior, which has implica- 
tions to a wide range of biological processes such as epige- 
netic inheritance, cell differentiation, and cancer oncoge- 
nesis • The nature and the existence of the multiple 
states are in the nonlinear biochemical reaction schemes 
and rates. For a traditional nonlinear bistable system in 
a small volume, the transition rates decreases with in- 
creasing system size. In the deterministic limit, these 
rates become zero. However, small systems can also ex- 
hibit bistable phenomenon which has no deterministic 
counter part. In the latter case, the bistability is due 
to the presence of "noise", "stochasticity" , or "entropic 
barrier". This class of bistable cellular systems can be 
quantitatively characterized by an van't Hoff-Arrhenius 
like analysis on the volume dependence of the transition 
rate(s). With increasing volume, the barrier diminishes. 
Mathematically, the transition rate is related to the land- 
scape of the stochastic dynamics, k oc e" v ^ x,v > where 
$(x, V) can be decomposed into <j>o(x) + <f>i/V. Nonlin- 
ear bistable system has barrier A^J > while stochastic 
bistability has A^j < but Ad>\ > 0. 
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